Compare AGB maps

Parameters

Year of ALOS-PALSAR mosaic: 2019 Backscatter processing variant: maskLWUwb_med5_LCinterp AGB variant: cap999pctl_maskWUwb

Look at distribution of AGB values

# Histogram of external map ----
bacc_hist <- plot_agb_hist_density(agb_fps$bacc, agb_pal)

## [1] "Sampling..."

hist_fp <- here::here('data/reports/external_agb', str_glue('agb_hist_Baccini.png'))
dir.create(dirname(hist_fp))
# ggsave(hist_fp, width = 5, height = 3, dpi = 120)

(esa_hist <- plot_agb_hist_density(agb_fps$esa, agb_pal))

## [1] "Sampling..."

hist_fp <- here::here('data/reports/external_agb', str_glue('agb_hist_ESA.png'))
# ggsave(hist_fp, width = 5, height = 3, dpi = 120)

(glob_hist <- plot_agb_hist_density(agb_fps$glob, agb_pal, xlim=NULL))

## [1] "Sampling..."

hist_fp <- here::here('data/reports/external_agb', str_glue('agb_hist_GlobB.png'))
# ggsave(hist_fp, width = 5, height = 3, dpi = 120)

(avit_hist <- plot_agb_hist_density(agb_fps$avit, agb_pal, xlim=NULL))

hist_fp <- here::here('data/reports/external_agb', str_glue('agb_hist_Avit.png'))
# ggsave(hist_fp, width = 5, height = 3, dpi = 120)

Compare overall coverage of each map

# Load data
smmry_csv <- here::here(comparison_dir, '07_overview_by_LC.csv')
stats_by_lc <- read_csv(smmry_csv) %>% 
  filter(!is.na(name))

# Helpers
lc_tbl <- list(`4` = 'Tree cover', `6` = 'Shrubs', `5` = 'Grassland', 
               `3` = 'Bareland', `2` = 'Urban', `1` = 'Water')
(map_order <- stats_by_lc %>% 
  mutate(map_code = factor(map_code, levels = c('internal', 'esa', 'glob', 'avit', 'bacc'))) %>% 
  arrange(desc(map_code)) %>% 
  distinct(name) %>% 
  deframe())

## [1] "Baccini"                          "Avitabile"                       
## [3] "GlobBiomass"                      "CCI"                             
## [5] "This study (cap999pctl_maskWUwb)"

# Prep DF
df_lc <- stats_by_lc %>% 
  filter(!is.na(name)) %>% 
  mutate(pct_of_lc_NA = 1 - pct_of_lc_100m) %>% 
  pivot_longer(cols = starts_with('pct_of_lc'),
               names_prefix = 'pct_of_lc_', 
               names_to = 'group', 
               values_to = 'pct') %>% 
  mutate(group = factor(group, levels = c('NA', '100m')), 
         name = factor(name, levels = map_order),
         LC = recode_factor(LC, !!!lc_tbl))

# Plot
(bp <- df_lc %>% 
    ggplot(aes(x = pct, y = name, fill = group)) +
    geom_col(position = 'fill', show.legend = FALSE) +
    scale_x_continuous(n.breaks = 3) +
    theme_minimal() +
    labs(fill="",x=NULL,y=NULL,title="",caption="") +
    facet_wrap(facets = vars(LC)))

# Save
ggsave(here::here(comparison_dir, '07_pct_masked_by_lc.png'),
       width = 6, height = 4, dpi = 120)

# Get percent of total
pct_total <- df_lc %>% 
    group_by(name, map_code) %>% 
    summarize(n_vals_100m = sum(n_vals_100m, na.rm = TRUE),
              n_lc_100m = sum(n_lc_100m, na.rm = TRUE)) %>% 
    dplyr::mutate(pct_of_lc_100m =  n_vals_100m / n_lc_100m,
           pct_of_lc_NA = 1 - pct_of_lc_100m) %>% 
  pivot_longer(cols = starts_with('pct_of_lc'),
               names_prefix = 'pct_of_lc_', 
               names_to = 'group', 
               values_to = 'pct') %>% 
  dplyr::mutate(group = factor(group, levels = c('NA', '100m')))

# Plot
(bp <- pct_total %>% 
    ggplot(aes(x = pct, y = name, fill = group)) +
    geom_col(position = 'fill', show.legend = FALSE) +
    scale_x_continuous(n.breaks = 3) +
    theme_minimal() +
    labs(fill="",x=NULL,y=NULL,title="",caption="") +
    ggtitle('Percent of all pixels (aggregated to 100m)'))

# Table
(d2 <- pct_total %>% 
    dplyr::mutate(pct_coverage = str_c(round(pct*100, digits = 1), '%')) %>% 
    dplyr::arrange(desc(pct)) %>% 
    dplyr::filter(!group %in% c('NA')) %>% 
    dplyr::select(name, pct_coverage))

tab <- gridExtra::tableGrob(d2, rows = NULL, cols = c('Map', 'Percent of land'))

# Display side-by-side
wrap_plots(tab, bp)

Compare external maps to our field data

# Scatterplot - AGB against backscatter ----------------------------------------
plot_scatter_ext <- function(plot_means, y_var, name){
  
  df <- plot_means %>% transmute(diff = .data[[y_var]] - AGB_ha) %>% deframe
  
  ggplot(plot_means, aes(x=AGB_ha, y=.data[[y_var]])) + 
    geom_point() +
    labs(x = expression(paste("Observed AGB (Mg ha"^"-1", ")")), 
         y = expression(paste("Estimated AGB (Mg ha"^"-1", ")")), 
         subtitle = str_c(name, ", mean difference: ", 
                          round(mean(df, na.rm=T), 1), " +/- ", 
                          round(sd(df, na.rm=T), 1), " t/ha")) +
    coord_cartesian(xlim=c(0,160), ylim=c(0,160)) +
    geom_abline(intercept = 0, slope = 1, linetype='dashed', alpha=.5) +
    geom_text(label=rownames(plot_means), hjust=0, vjust=0, size=3.5) + 
    theme_minimal()
  
}

# Load data
plot_means <- read_csv(plot_ext_csv)

# Extract just AGB and backscatter as dataframe
p0 <- plot_scatter_ext(plot_means, 'internal', name='Our AGB')
p1 <- plot_scatter_ext(plot_means, 'glob', name='GlobBiomass')
p2 <- plot_scatter_ext(plot_means, 'esa', name='ESA')
p3 <- plot_scatter_ext(plot_means, 'avit', name='Avitabile')
p4 <- plot_scatter_ext(plot_means, 'bacc', name='Baccini')

# Save
plot_fp <- here::here(comparison_dir,
                     str_c('comparison_scatter_', agb_code, '_ourAGB.png'))
ggsave(plot_fp, p0, width=10, height=10, units='cm')

(p1 | p2 )/ (p3 | p4)

plot_fp <- here::here(comparison_dir,
                     str_c('comparison_scatter_', agb_code, '.png'))
ggsave(plot_fp, width=20, height=20, units='cm')

Compare pixel-to-pixel difference metrics (R^{2}, RMSD, Bias)

# Load
sum_tbl <- read_csv(compare_by_lc_csv) %>% 
  filter(!is.na(R2), 
         mask != 'WUBGS')

# Experiment with views
long_tbl <- sum_tbl %>% 
  rename(product = name) %>% 
  pivot_longer(R2:med_diff, names_repair = 'unique')

long_tbl %>% 
  filter(name %in% c('R2', 'RMSD', 'bias')) %>% 
  ggplot(aes(x = product, y = value, color = product)) +
  geom_point() +
  facet_grid(rows = vars(name), cols = vars(mask), scales = 'free_y') +
  theme(axis.title = element_blank())

long_tbl %>% 
  filter(name %in% c('MBD')) %>% 
  ggplot(aes(x = product, y = value, color = product)) +
  geom_point() +
  geom_hline(yintercept = 0) +
  facet_grid(rows = vars(name), cols = vars(mask), scales = 'free_y') +
  theme(axis.title = element_blank())

long_tbl %>% 
  filter(name %in% c('R2', 'RMSD', 'bias'), 
         mask %in% c('TC', 'BGS')) %>% 
  ggplot(aes(x = mask, y = value, color = product)) +
  geom_point() +
  facet_wrap(vars(name), scales = 'free_y') +
  theme(axis.title = element_blank())

Compare pixel-to-pixel scatter density plots

# Tree cover
mskvals <- mskvals_list[[1]]
agb_fps_sub <- agb_fps[c(2,3,4,5)]
  
plots_tc <- agb_fps_sub %>% 
      purrr::map(compare_by_given_lc, agb_fp, lc_fp,
                  mskvals = mskvals, comparison_dir,
                  return_obj = 'plot', 
                  boundary_fp = hti_poly_fp) %>% 
    flatten()

plots_tc %>% saveRDS(here::here(comparison_dir, 'plot_tc.rds'))

# Plots for given mask
(patches <- plots_tc %>% 
  wrap_plots(nrow = 4, guides = 'collect') &
    theme(axis.title = element_blank()))
map_fp <- here::here(comparison_dir, str_glue("scattdens_{mskvals$code}.png"))
ggsave(map_fp, width=4, height=12, dpi=120)

# Total 
mskvals <- mskvals_list[[2]]
agb_fps_sub <- agb_fps[c(2,3,4,5)]
  
plots_total <- agb_fps_sub %>% 
      purrr::map(compare_by_given_lc, agb_fp, lc_fp,
                  mskvals = mskvals, comparison_dir,
                  return_obj = 'plot', 
                  boundary_fp = hti_poly_fp) %>% 
    flatten()

plots_total %>% saveRDS(here::here(comparison_dir, 'plots_total.rds'))

# Plots for given mask
(patches <- plots_total %>% 
  wrap_plots(nrow = 4, guides = 'collect') &
    theme(axis.title = element_blank()))
map_fp <- here::here(comparison_dir, str_glue("scattdens_{mskvals$code}.png"))
ggsave(map_fp, width=4, height=12, dpi=120)

# BGS 
mskvals <- mskvals_list[[3]]
agb_fps_sub <- agb_fps[c(2,3,4,5)]
  
plots_bgs <- agb_fps_sub %>% 
      purrr::map(compare_by_given_lc, agb_fp, lc_fp,
                  mskvals = mskvals, comparison_dir,
                  return_obj = 'plot', 
                  boundary_fp = hti_poly_fp) %>% 
    flatten()

plots_bgs %>% saveRDS(here::here(comparison_dir, 'plots_bgs.rds'))

# Plots for given mask
(patches <- plots_bgs %>% 
  wrap_plots(nrow = 4) &
    theme(axis.title = element_blank()))
map_fp <- here::here(comparison_dir, str_glue("scattdens_{mskvals$code}.png"))
ggsave(map_fp, width=4, height=12, dpi=120)

# Combine

(patches <- c(plots_tc, plots_bgs, plots_total) %>% 
  wrap_plots(nrow = 4, ncol = 3, 
             byrow = FALSE,
             guides = 'collect') &
    theme(axis.title = element_blank()))

map_fp <- here::here(comparison_dir, str_glue("scattdens_all.png"))
ggsave(map_fp, width=12, height=12, dpi=120)


plots_tc$glob +
  guides(color = guide_legend('density'),
         alpha = guide_legend('density'))

plots_tc$glob + guides(alpha = FALSE)

# Tree cover
mskvals <- mskvals_list[[1]]
agb_fps_sub <- agb_fps[c(2,3,4,5)]
  
df <- compare_by_given_lc(agb_fps_sub[[1]], agb_fp, lc_fp,
                  mskvals = mskvals, comparison_dir,
                  return_obj = 'diffs', 
                  boundary_fp = hti_poly_fp)
ext_name = ''
sample_size = 50000
xlim = c(0, 300)



# Guides are integrated where possible
p + guides(colour = guide_legend("title"), size = guide_legend("title"),
  shape = guide_legend("title"))


# ggplot object
dat <- data.frame(x = 1:5, y = 1:5, p = 1:5, q = factor(1:5),
 r = factor(1:5))
p <- ggplot(dat, aes(x, y, colour = p, size = q, shape = r)) + geom_point()

# without guide specification
p

# Show colorbar guide for colour.
# All these examples below have a same effect.
p + guides(colour = "colorbar", size = "legend", shape = "legend")
p + guides(colour = guide_colorbar(), size = guide_legend(),
  shape = guide_legend())
p +
 scale_colour_continuous(guide = "colorbar") +
 scale_size_discrete(guide = "legend") +
 scale_shape(guide = "legend")

 # Remove some guides
 p + guides(colour = "none")
 p + guides(colour = "colorbar",size = "none")

# Guides are integrated where possible
p + guides(colour = guide_legend("title"), size = guide_legend("title"),
  shape = guide_legend("title"))

# same as
g <- guide_legend("title")
p + guides(colour = g, size = g, shape = g)

p + theme(legend.position = "bottom")

# position of guides

# Set order for multiple guides
ggplot(mpg, aes(displ, cty)) +
  geom_point(aes(size = hwy, colour = cyl, shape = drv)) +
  guides(
   colour = guide_colourbar(order = 1),
   shape = guide_legend(order = 2),
   size = guide_legend(order = 3)
 )

Look at AGB maps

plot_agb_map <- function(x, mask_poly = NULL, borders = NULL, resolution = 200) {
  
  lyr_name <- x$name
  
  # Load AGB
  agb <- rast(x$fp)
  
  # Conditionally aggregate map based on resolution
  resx <- res(agb) %>% nth(1) #%>% signif(2)
  fact <- resolution/resx/113000
  agg_factor <- ifelse(fact >= 1.5, round(fact, 0), 1)
  if (agg_factor > 1) {
    agb <- terra::aggregate(agb, fact = agg_factor, fun = 'mean', na.rm = TRUE)
  } else if(fact < 1) {
    print('WARNING: input map resolution is greater than target')
  }

  # Convert raster to dataframe
  agb_a_df <- agb %>% as.data.frame(xy = T) %>% rename(agb = 3)

  if(is.null(mask_poly)){
      # Urban polygons
      tol <- 0.001
      urb_poly_fp <- here::here(tidy_dir, 'landcover', str_glue('lc2_urban_tol{tol}.gpkg'))
      if(!file.exists(urb_poly_fp)) {
        mask_poly <- st_read(lc_pols_fp) %>% 
          filter(LC == 2) %>% 
          st_simplify(dTolerance = tol) %>% 
          st_union()
        mask_poly %>% st_write(urb_poly_fp)
      } else {
        mask_poly <- st_read(urb_poly_fp)
      }
  }

  if(is.null(mask_poly)){
    # Haiti administrative boundaries
    borders <- st_read(here::here(tidy_dir, 'contextual_data', 'HTI_adm', 'HTI_adm1.shp'))
  }
  
  # Plot
  agb_map <- ggplot(agb_a_df) +
    geom_raster(aes(x = x, y = y, fill = agb)) +
    # geom_tile(aes(x = x, y = y, fill = agb)) +
    geom_sf(data = mask_poly, fill = "gray63", color = NA) +
    geom_sf(data = borders, fill = NA, color = "gray39", size = 0.3) +
    scale_fill_gradientn(expression(paste("Aboveground\nbiomass (Mg ha"^"-1", ")")), 
                         na.value = "transparent", 
                         colors = agb1_palette,
                         limits = c(0, 120),
                         oob = scales::squish
    ) +
    ggthemes::theme_map() +
    ggtitle(lyr_name) +
    theme(legend.position=c(0.23, 0.9),
          legend.title.align=0,
          legend.justification = c(1,1),
          legend.box.background = element_blank(),
          panel.grid.minor = element_blank(),
          axis.title = element_blank(), 
          plot.title = element_text(vjust = -8, hjust = 0.048)) 
}

# Urban polygons
tol <- 0.001
urb_poly_fp <- here::here(tidy_dir, 'landcover', str_glue('lc2_urban_tol{tol}.gpkg'))
if(!file.exists(urb_poly_fp)) {
  urban <- st_read(lc_pols_fp) %>% 
    filter(LC == 2) %>% 
    st_simplify(dTolerance = tol) %>% 
    st_union()
  urban %>% st_write(urb_poly_fp)
} else {
  urban <- st_read(urb_poly_fp)
}

## Reading layer `lc2_urban_tol0.001' from data source `/home/esturdivant/PROJECTS/biomass-espanola/data/tidy/landcover/lc2_urban_tol0.001.gpkg' using driver `GPKG'
## Simple feature collection with 1 feature and 0 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.4575 ymin: 18.04225 xmax: -71.68825 ymax: 20.04025
## Geodetic CRS:  WGS 84

# Haiti administrative boundaries
hti_adm1 <- st_read(here::here(tidy_dir, 'contextual_data', 'HTI_adm', 'HTI_adm1.shp'))

## Reading layer `HTI_adm1' from data source `/home/esturdivant/PROJECTS/biomass-espanola/data/tidy/contextual_data/HTI_adm/HTI_adm1.shp' using driver `ESRI Shapefile'
## Simple feature collection with 10 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.48125 ymin: 18.0218 xmax: -71.61815 ymax: 20.09042
## Geodetic CRS:  WGS 84

Our AGB map

resolution <- 400

(agb_map <- plot_agb_map(agb_fps$internal, urban, hti_adm1, resolution) +
   theme(plot.title = element_blank()))

Data resampled to 400 m resolution for plotting.

This study, GlobBiomass, ESA, and Baccini

# Extract plot means for all AGB maps
agb_fps$internal$name <- 'This study'
maps <- agb_fps %>% 
  purrr::map(plot_agb_map, urban, hti_adm1, resolution = resolution) 

## [1] "WARNING: input map resolution is greater than target"

# Set output dir for side-by-side maps
out_dir <- here::here(comparison_dir, 'sidebyside')
dir.create(out_dir)

# Ours, with GlobBiomass, ESA, and Baccini
all_maps <- maps[c(1,3,2,5)] %>% 
  wrap_plots(nrow = 2, guides = 'collect')
map_fp <- here::here(out_dir, str_glue("sidebyside_Us_ESA_GlobB_Bacc_{resolution}m.png"))
ggsave(map_fp, width=12, height=10, dpi=120)

all_maps

Avitabile, ESA, GlobBiomass, and Baccini

# 4 external maps
all_maps <- maps[c(4,2,3,5)] %>% 
  wrap_plots(nrow = 2, guides = 'collect')
map_fp <- here::here(out_dir, str_glue("sidebyside_Avit_GlobB_ESA_Bacc_{resolution}m.png"))
ggsave(map_fp, width=12, height=10, dpi=120)

all_maps

This study, GlobBiomass, Avitabile, ESA, and Baccini

# All 5, equal sizes
maps$space <- guide_area()
all_maps <- maps[c(1, 3, 2, 6, 4, 5)] %>% 
  wrap_plots(guides = 'collect')
map_fp <- here::here(out_dir, str_glue("all5_equal_{resolution}m.png"))
ggsave(map_fp, width=14, height=9, dpi=120)

all_maps

This study, GlobBiomass, Avitabile, ESA, and Baccini

# All 5, with ours bigger
four <- ((maps$bacc + maps$avit) / (maps$glob + maps$esa)) &
  theme(legend.position = 'none')
all_maps <- (maps$internal + plot_layout(guides = 'keep')) / four
map_fp <- here::here(out_dir, str_glue("all5_{resolution}m.png"))
ggsave(map_fp, all_maps, width=5, height=8, dpi=120)

all_maps

Difference maps

# ~ Difference maps (only for all LCs) -----------------------------------------
mskvals <- list(values = c(1,2,3,4,5,6), code = 'Total')
diffmaps <- agb_fps[c(2,3,4,5)] %>% 
  purrr::map(compare_by_given_lc, agb_fp, lc_fp,
              mskvals = mskvals, comparison_dir,
              return_obj = 'map', 
              boundary_fp = hti_poly_fp)

# 4 external maps
all_maps <- diffmaps %>% 
  wrap_plots(nrow = 2, guides = 'collect') &
    ggthemes::theme_map()
map_fp <- here::here(out_dir, str_glue("difference_maps.png"))
ggsave(map_fp, width=14, height=10, dpi=120)

all_maps